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We study spherical and nonspherical linear acoustic perturbations of the Michel flow, which 
describes the steady radial accretion of a perfect fluid into a nonrotating black hole. The dynamics 
of such perturbations are governed by a scalar wave equation on an effective curved background 
geometry determined by the acoustic metric, which is constructed from the spacetime metric and 
the particle density and four-velocity of the fluid. For the problem under consideration in this article 
the acoustic metric has the same qualitative features as an asymptotically flat, static and spherically 
symmetric black hole, and thus it represents a natural astrophysical analogue black hole. 

As for the case of a scalar field propagating on a Schwarzschild background, we show that acoustic 
perturbations of the Michel flow exhibit quasi-normal oscillations. Based on a new numerical method 
for determining the solutions of the radial mode equation, we compute the associated frequencies 
and analyze their dependency on the mass of the black hole, the radius of the sonic horizon and 
the angular momentum number. Our results for the fundamental frequencies are compared to those 
obtained from an independent numerical Cauchy evolution, finding good agreement between the two 
approaches. When the radius of the sonic horizon is large compared to the event horizon radius, we 
find that the quasi-normal frequencies scale approximately like the surface gravity associated with 
the sonic horizon. 

PACS numbers: 04.20.-q,04.70.-s, 98.62.Mw 


I. INTRODUCTION 

The study of accretion into a black hole plays a very important role in general relativity and astrophysics. In 
particular, an understanding of the emission of electromagnetic radiation generated by compression or friction in the 
gas is an important subject since this radiation may carry information about the spacetime geometry close to the 
black hole and thus offer the opportunity to test Einstein’s general theory of gravity in its strong field limit. In fact, 
millimeter-wave very-long baseline interferometric arrays such as the Event Horizon Telescope flj are already able 
to resolve the region around Sagittarius A*, the supermassive black hole lying in the center of our galaxy, to scales 
smaller than its gravitational radius [2], Comparing the observations to calculated images of the black hole shadow 
and the sharp photon ring surrounding it may even lead to tests for the validity of the no-hair theorems [3]. 

Clearly, the features of the observed electromagnetic signals depend on the properties and dynamics of the flow, 
and therefore it is of considerable interest to study the dynamics of the accreted gas and to identify its key properties 
like its oscillation modes, for example. For a numerical study of oscillating relativistic fluid tori around a Kerr black 
hole and astrophysical implications, see Ref. j4j . For the impact of a binary black hole merger on the dynamics of the 
circumbinary disk and associated electromagnetic signals, see Refs. mm and references therein. 

Motivated by the above considerations, the purpose of the present work is to study the oscillation modes of a 
simple accretion model, namely the radial flow of a perfect fluid on a nonrotating black hole background. Spherically 
symmetric steady-state configurations in this model for which the density is nonzero and the matter is at rest at 
infinity have been studied long time ago by Michel 0, generalizing previous work by Bondi 0 in the Newtonian case. 
The Michel flow describes a transonic flow, the flow’s radial velocity measured by static observers being subsonic in 
the asymptotic region and supersonic close to the event horizon. Although much less realistic than the case where 
the black hole rotates and/or the matter has an intrinsic angular momentum, resulting in an accretion disk, the 
study of spherical accretion is still relevant in a variety of interesting astrophysical scenarios. Examples include 
nonrotating black holes accreting matter from the interstellar medium 01 and supermassive black holes accreting 
dark matter [10]. For a rigorous treatment on the Michel flow and its generalization to a wide class of spherical black 
hole backgrounds, we refer the reader to our recent work m- 

In this article, we study spherical and nonspherical linear acoustic perturbations of the Michel flow, assuming a 
fixed Schwarzschild black hole background. Moncrief [T2j showed that if the entropy and vorticity perturbations are 
of bounded extent on some initial hypersurface, they will be advected into the black hole in finite time, leaving a 
pure potential flow perturbation in their wake. Furthermore, Moncrief showed in Ref. [ T21 that the potential flow 
perturbation can be described in a very elegant manner by a wave equation on an effective curved background geometry 
described by the acoustic (or sound) metric , which is constructed from the spacetime metric and the four-velocity 
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and particle density of the background flow. The acoustic metric is Lorentzian and its null cones (the sound cones) 
lie inside the light cones, as long as the speed of sound is smaller than the speed of light. For further properties of 
the acoustic metric, see [151 . 

For acoustic perturbations of the Michel flow the geometry described by the acoustic metric is asymptotically flat, 
static and spherically symmetric and possesses a sonic horizon, defined as the boundary of the region which can send 
sound signals to a distant observer, where the matter is almost at rest. As it turns out, this boundary coincides with 
the location of the sonic sphere describing the transition of the flow’s radial velocity measured by static observers from 
sub- to supersonic. Therefore, as far as the propagation of sound waves are concerned, the acoustic geometry for the 
Michel flow has exactly the same qualitative properties as the geometry of a static, spherically symmetric black hole 
on which electromagnetic radiation propagates, and the sonic horizon in the acoustic geometry plays the role of the 
event horizon. Consequently, the acoustic geometry for the Michel flow constitutes a natural astrophysical “analogue 
black hole”. For a review on analogue black holes in different physical contexts, we refer the reader to Ref. m, and 
for recent applications to accretion flows on black hole backgrounds, see Refs. [TSHU]. 

Interpreting the acoustic perturbations as an evolution problem on an effective geometry leads to new insight and 
new results. For the case of the Michel flow, for example, one can prove that acoustic perturbations outside the 
sonic horizon stay bounded, using standard energy conservation techniques nmnmnBi- In this article, we use this 
analogue black hole interpretation and show that, similar to the case in which a Schwarzschild black hole is perturbed, 
small perturbations of the Michel flow lead to quasi-normal acoustic oscillations characterized by complex frequencies 
s = a + ioj, where a < 0 describes the decay rate and ui the frequency of oscillation. As in the black hole case, 
these frequencies describe the ringdown phase which is taken over by a slower power-law decay at late times. We 
numerically compute the quasi-normal frequencies (and in some cases also the exponent in the late-time power-law 
tail) as a function of the black hole mass (or its Schwarzschild radius ?'#), the radius of the sonic horizon r c and the 
angular momentum number £ of the perturbation. For previous studies of quasi-normal oscillations in fluid analogue 
black hole modes, see for example Refs. mum no- Contrary to these references which are mainly concerned with 
analogue black holes in the laboratory, the scenario considered in this article refers to an astrophysical analogue black 
hole. 

The remainder of this work is organized as follows. In Sec. |TT] we briefly review the main features of the Michel 
flow, and in particular we discuss the properties of the flow in the vicinity of the sonic sphere. Next, in Sec. Ill 


we first analyze the geometric properties of the acoustic metric and show that it indeed describes an analogue black 
hole whose horizon is located at the sonic sphere. We also compute the surface gravity associated with this sonic 
horizon since it plays an important role in the description of the quasi-normal acoustic frequencies found in this article. 
Next, by performing a mode decomposition, we reduce the wave equation on the acoustic metric background to a 
family of radial, time-independent Schrodinger-like equations and discuss our method for computing the quasi-normal 
frequencies. One important issue we would like to point out here is that unlike the case where the background metric 
is Schwarzschild, the effective potential appearing in our radial equation cannot be written in explicit form. This 
complication stems from the fact that the Michel solution, describing the particle density as a function of the areal 
radius coordinate, is only known in implicit form, and consequently the metric coefficients in the acoustic metric and 
the effective potential in the radial equation can only be described in terms of implicit functions. For this reason the 
problem is much harder than in the Schwarzschild case, and popular analytic methods based on series expansions 
like Leaver’s method [22] do not seem feasible. This issue has motivated us to reconsider the problem of calculating 
the quasi-normal frequencies based on a new numerical matching procedure, where the local solutions of the radial 
equation which are being matched are computed via a Banach iteration method. This method, which shares some 
common features with the complex coordinate WKB approximation (see [55] and references therein), is described and 
tested in Sec. |III| See also [24] for a recent method allowing to compute the quasi-normal frequencies for deformed 


Kerr black holes based on ideas from perturbation theory in quantum mechanics. 

Next, in Sec. |IV| we describe a completely different method for computing the quasi-normal frequencies based on a 
numerical Cauchy evolution of the wave equation. In this method, one specifies an initial perturbation for the fluid’s 
acoustic potential, solves the wave equation numerically and registers the signal observed by a static observer outside 
the sonic horizon. The signal reveals an initial burst followed by a ringdown signal whose oscillations frequency uj 
and decay rate a can be combined into a complex frequency. Comparing s = a + iu with the fundamental quasi¬ 
normal frequencies computed in Sec. |III| provides a further validation for our matching procedure, and shows that the 
quasi-nornral acoustic oscillations found in this paper are actually excited by an initial perturbation of the fluid. The 
numerical results in Sec. |IV| also indicate that the ringdown signal is overtaken by a power-law decay at late times, 
similar to what has been observed in laboratory-type analogue black holes pH] . 

Our main results for the quasi-normal acoustic frequencies and their dependency on r c and th and on the angular 
momentum number t are presented in Sec. [V]for the case of a polytropic fluid equation of state with adiabatic index 
7 = 4/3. Our results indicate that for large r c /r# the complex frequencies s scale approximately like the surface 
gravity n of the acoustic geometry and that the rescaled decay rates a/n do not depend strongly on £ for r c r# 
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and £ > 1. Results for overtone frequencies and the eikonal limit i —> oo are also discussed in Sec. [Vj Conclusions are 
drawn in Sec. |VI| and technical details related to the analytic continuation of the effective potential needed for our 
matching procedure are explained in an appendix. 


II. REVIEW OF MICHEL FLOW AND ITS RELEVANT PROPERTIES 

In this section, we review the relevant equations describing the Michel flow on a Schwarzschild background. For 
details and a generalization to more general static, spherically symmetric black hole backgrounds, see Refs. Q31I251 [26]. 
We write the Schwarzschild metric in the form 

g = — N{r)c 2 dt 2 + fj" r 2 (dfl 2 _|_ s j n 2 fldip 2 ) ; N(r) = 1 — —, (1) 

N(r) r 

where c is the speed of light and th the Schwarzschild radius. The fluid is described by the particle density n, energy 
density e and pressure p measured by an observer moving along the fluid four-velocity u = u^d^. (u is normalized 
such that = — c 2 .) Its dynamics is determined by the equations of motion 


V M = 0, 


( 2 ) 


= 0 . 


( 3 ) 


wherein = nu M is the particle current density and T= nhiFiF + pg IJjl/ is the stress-energy tensor, and V refers 
to the covariant derivative with respect to the spacetime metric g. Here and in the following, h denotes the enthalpy 
per particle, defined as h := (p + s)/n , and we assu me that h = h(n) is a function of the particle density n only. In 
the spherically symmetric stationary case Eqs. ( |2|3[ ) reduce to 


47 rr 2 ?ru r = j n = const, 


4 - 7 T r 2 nhu r \ N + 


= j e = const, 


( 4 ) 

( 5 ) 


which expresses the conservation of particle and energy flux through a sphere of constant areal radius r. Using Eq. Q 
in order to eliminate u r in Eq. |5]) gives 


F(r, n ) 


h{n) 


N(r ) 


,2 1 



= const, 


jn p. 

' ,:= te <0 ' 


( 6 ) 


where j n describes the accretion rate and is negative. Therefore, the problem of determining the accretion flow is 
reduced to finding an appropriate level curve of the function F(r,n ), which associates to each value of r a unique 
value of the particle density n(r). Once n(r) is known, the radial velocity u r is obtained from Eq. Q. 

In previous work m we proved that under the conditions on the equation of state (F1)-(F3) below there exists 
a unique smooth solution n(r ) of Eq. Q which extends from the event horizon r = rjj to infinity and has a given 
positive particle density rioo > 0 at infinity. We shall call this solution the Michel solution. Our conditions on h(n), 
which we assume to be a smooth function h : (0, oo) —> (0, oo), are the following: 

(FI) h{n) — > e-o > 0 for n -A 0 (positive rest energy), 

(F2) 0 < = q < 1 for all n > 0 (positive and subluminal sound velocity), 

(F3) 0 < W(n) := | for all n > 0 (technical restriction on the derivative of v s (n)). 

In particular, these conditions are satisfied for a polytropic equation of state 

h{n) = e 0 + Kn y ~ 1 , (7) 


wherein eo > 0, K > 0 and the adiabatic index 7 lies in the range 1 < 7 < 5/3. In this work, we focus on the 
particular case of an ultra-relativistic gas for which h(n) has the same form as in Eq. 0 with 7 = 4/3. However, for 
the sake of generality, all the expressions below are given for an arbitrary equation of state satisfying the assumptions 
(F1)-(F3). 
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The function n : [rn, oo) —>• R. describing the Michel flow is a smooth, monotonously decreasing function which is 
implicitly determined by Eq. ([6]), that is 

F(r,n(r)) = const = h(n oo) 2 > 0. 


By differentiating both sides with respect to r one obtains 


r,n(r )) + ^-(r,n(r))n'(r) = 0, 


where the partial derivatives of F are 


dF . h(n) 2 rn 4/r 2 

— (r,n = - 

or r 


dF , i 

aU’ n) = 


2 h(n) 2 


2 r .,2 


-J-Af(r) — ( 1 — -4 


,2 1 


( 8 ) 

( 9 ) 

( 10 ) 


The implicit function theorem guarantees local existence and uniqueness of n(r) as long as dF/dn ^ 0. In the 
asymptotic region (large r), dF/dn > 0 is positive, and close to the event horizon (r ~ m) dF/dn < 0 is negative, 
so in these regions the slope n' of n is uniquely determined by Eq. ([8]). However, by continuity, there exists a point 
r c > th where dF/dn vanishes, and at this point n'[r c ) can only be finite if dF/dr also vanishes. This leads to the 
requirement that the flow must necessarily pass through a critical point (r c ,n c ) of the function F{r,n). In [TTi we 
proved that under the assumptions (F1),(F2),(F3) on the fluid there is, for large enough |/z|, a unique critical point of 
F(r, n) and a unique solution n(r) of Eq. (Jbj) which extends from ru to r = oo and satisfies n(r c ) = n c . Furthermore, 
given rioo > 0 the value of |/x| (and hence the location of the critical point) is fixed. 

Physically, the critical point corresponds to the sonic sphere r = r c , which describes the transition of the flow’s 
radial velocity measured by static observers from sub- to supersonic. The location of the sonic sphere is determined 
by the equations 


— = t I 3+ — 
r H 


M 


V s 

V c ■= ~ , 
c 


( 11 ) 


which follow from setting the right-hand sides of Eqs. ( 9|l0 ) to zero. According to assumption (F2), v ~ 2 is always 
larger than one, and Eq. (11) implies that the sonic horizon is located outside the horizon. 

For later use we shall also need the derivative of the particle density n' c := n'{r c ) at the critical point. For this, we 
differentiate Eq. (18)) with respect to r and evaluate at r = r c , obtaining 


d 2 F, , „ d 2 F , . , d 2 F 

n c ) + 2 ^-( r c, n c )n c + — (r c 


dr 2 ’ drdn 

Using the following expression for the Hessian of F at (r c ,n c ), 

K 2 


dn 2 


o)K) 2 = o. 


( 12 ) 


gr -2 (v c ,n c ) g r f) n (j~C1 n c) 
££(r c ,n c ) 0(r c ,n c )y 

with h c = h{r c ) and W c = W(r c ), we find the two solutions 



Vc 

n c 


r c 2±^l + 3{v 2 ~W c y 


(13) 


which parametrize the two branches of the level set of F through (r c ,n c ). In m we proved that the branch corre¬ 
sponding to the global solution for n(r) extending from the horizon to infinity is the one with the + sign in Eq. (13). 
In the appendix, we show that the function n(r) admits an analytic continuation to complex r. This continuation is 
required for the quasi-normal mode calculation in the next section. 


III. QUASI-NORMAL OSCILLATIONS FROM A MODE ANALYSIS 


The propagation of acoustic perturbations in any relativistic perfect fluid is elegantly described by a wave equation 

□e* = 0, (14) 
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where the scalar field U/ determines the perturbed enthalpy Sh and four-velocity Su M of the fluid according to the 
relation S(hu M ) = V^, which using u^Su^ = 0 yields 

8h = -u fi V^, ^ = 

The operator □© in Eq. ( |14[ ) is the wave operator belonging to the acoustic metric 0, which is constructed from the 
spacetime metric g and the fluid quantities in the following way m- 


<v 


n c 
h v s 


9u" 



(15) 


Under our assumptions on the sound speed it follows that 0 is a Lorentzian metric whose cone (the sound cone) lies 
inside the light cone of g. Notice also that u is timelike with respect to both g and 0. 


A. Geometry of the acoustic metric 


For simplicity, from now on we use units in which the speed of light is one, c = 1. For the particular case of the 
Michel flow on a Schwarzschild metric the acoustic metric is 


© = 


n 1 
h v s 


dr 2 

—Ndt 2 + -^r- + (1 — v 2 ) (Utdt + u r dr ) 2 + r 2 (dr? 2 + sin 2 ddip 2 ) 


or 


with 


0 = — A(r)dt 2 + 2 B{r)dtdr + C(r)dr 2 + R(r) 2 (dr? 2 + sin 2 r?dyj 2 ) , 


(16) 


(17) 


A{r) 
B(r ) 
C(r) 
R(r) 


2l[„ 2A r_ ( i 

n Vo 


n 1 — v 2 y/N + ( u r ) 2 u r 
h v s 


N 


n 1 1 
hmN 2 


[iV + (1 — v 2 )(u r ) 2 


n 1 

1 r ; 
h Vs 


where we have used the equation u 2 — ( u r ) 2 = N and Ut < 0 in order to eliminate ut and where the quantities n, h, v s 
and u r are given by the Michel flow solution discussed in the previous section. The acoustic metric (16) is spherically 
symmetric and possesses the Killing vector field 


k -l 

dt 


(18) 


whose negative square norm is A{r). Since A(r) is positive for r > r c and negative for r < r c (cf. Eq. (101 and 
the remarks following this equation) the vector field k is timelike for r > r c , spacelike for 0 < r < r c and null at 
r = r c , and the surface r = r c is a Killing horizon HOT- Notice that the coordinates (f, r) are regular everywhere 
outside the event horizon r > r#; in particular they are regular at the sonic horizon r = r c . Introducing the new time 
coordinate 


T \=t — 


B (r) 

A{r) 


dr , 


the acoustic metric can be brought into diagonal form outside the sonic horizon, 


0 = 


n 1 

h v s 


dr 


-X{r)v 2 dT 2 + + r 2 (di? 2 + sin 2 r?d^ 2 ) 


X(r) := A(r)- --1 (rU) 2 . 


(19) 


Note that X(r) —► 1 as r —> oo, and in this limit the acoustic metric reduces (up to a constant conformal factor) to 
the Minkowksi metric with time coordinate Coo T, where Vac := linv^oo v s (r) is the sound speed at infinity. 
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It follows that the geometry described by the acoustic metric (16) is the same as the one of a static, spherically 
symmetric and asymptotically flat black hole. The sonic horizon r = r c plays the role of the event horizon of this 
analogue black hole. Its surface gravity n with respect to the Killing vector field k defined in Eq. (18), which will play 
an important role later, can be computed using Eqs. (Ill and (131. The result is 


^ A'(r c ) 
2 B(r c ) 


1 r H 
Av c r 2 


Vi + W! 


w c ). 


( 20 ) 


Since scalar fields propagating on static spherically symmetric black holes like the Schwarzschild and Reissner- 
Nordstrom black holes exhibit quasi-normal oscillations, and since the fluid potential T satisfies a wave equation on 
an analogue black hole background, it is natural to expect that acoustic perturbations in the Michel flow undergo 
quasi-nornral oscillations as well. In the following, we show that such oscillations do indeed exist and compute the 
associated frequencies based on two different numerical methods. 


B. Reduction to a Schrodinger-like equation 


Quasi-normal modes are particular solutions of Eq. (14) which are of the form 


sT 1>{s,r)Y tm {-d,<p\ 

for some complex frequency s = a + iu> € C and complex-valued function ip(s,r) to be determined. Here, a denotes 
the decay rate, u> the frequency of oscillations, and Y fm the standard spherical harmonics with angular momentum 
numbers £m. Introducing this ansatz into Eq. (14) and using the diagonal parametrization (191 of the acoustic metric, 
one obtains the following equation: 




\r, 


+ [s 2 + Af(r)Ve(r)\ ip = 0, 


( 21 ) 


where the functions A f(r) and Vi(r) are explicitly given by 


A f(r) = v s X = v , 

1 


1 - 


rn 


-iK ) 2 


( 22 ) 


Vt(r) = 


r n 


-(1 - V i + 5W)E-~- 


2 n 


AW + 3W 2 + (1 - vi) 2 - 2 


dW 

d\ogn_ 


E+ K 
4r 


1 + 3vi + 3W + AWr— 
s n 


K + i) 


(23) 


with E := rn/{Ar) — ( u r ) 2 , u r = [i/(r 2 n ) and W = d\ogv s /dlogn. Away from the critical point n'/n can be 
computed using Eq. ([8]), which yields 


n 

n 


2 E 


r v 2 X ’ 


(24) 


while for r = r c Eq. (13) can be used in order to compute n'/n. 

For large r the effective potential Vi{r) behaves as u 00 £(£ + l)/r 2 + 0(r -3 ), so it is dominated by the centrifugal 
term. At the sonic horizon A f(r c ) is zero, but V)(r c ) is positive. Introducing the tortoise coordinate r» = f dr/Af(r) 
which ranges from — oo to +00 Eq. (21) can be further simplified and is formally equivalent to the time-independent 
Schrodinger equation Hip = — s 2, ip with Hamiltonian 


H := -^2 +N{r)Vt(r). 

A plot of the effective potential J\f(r)Ve(r) for t = 0 is shown in Fig. [T] which indicates that there is a potential barrier 
even in the monopolar case £ = 0. 
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W(x) 



FIG. 1: The effective dimensionless potential W(x) := r| f A/”(r)Vb(r) in the Hamiltonian H as a function of x := t/xh is shown 
here for the case of the Michel flow for a polytrope with adiabatic index 7 = 4/3 and sonic horizon located at r c = 2xh- 


C. Computation of the quasi-normal frequencies using Banach iterations 


The quasi-normal frequencies s are determined by the following requirement [29j . 
admits precisely two solutions 4’±( s i r ) satisfying the boundary conditions 

lim e sr *ip+( s , r) = 1 , lim e~ sr *ip-(s, r) = 1 , 

r*—>oo r *—>■—00 


For a = Re(s) > 0 Eq. (21) 


(25) 


in the asymptotic region and at the sonic horizon, respectively. These solutions can be shown to depend analytically 
on s, and they can be analytically continued on the left complex plane a < 0. For a > 0 the two functions ip+(s, •) 
and ip-( s j ') are always linearly independent from each other since otherwise one would have a finite energy solution 
which grows exponentially in time, in contradiction to standard energy arguments | 12 | showing the stability of the 
flow outside the sonic horizon. However, for particular values of the complex frequency s = a + ioj with a < 0 it is 
possible that the two functions become linearly dependent. These special frequencies are the ones associated with 
the quasi-normal modes, and as we will show in the next section they describe the ringdown phase in the dynamics 
of the scalar field ’F. For more general discussions on quasi-normal oscillations we refer the reader to the review 
articles [5DH52] , 

For Re(s) > 0 the solutions if>± can be constructed using the following iteration scheme, 

i>±(s,r) = e Tsr * lim (T^ s l)(r), (26) 

k—> 00 

where the operators T± s , acting on continuous and bounded functions £, are defined as 

OO 

(T +a Q(r) = 1 +^| (l - e - 2 «(rl-r.)) Vt (r , )^r')dr f , 

r 

r 

(T-e$){r) = 1+^/ (l - e+ 2 ^*-*•*>) V e (r'W)dr', 

r c 

for r c < r < 00, with r' the variable of integration and r* the associated tortoise coordinate. Note that the integrals 
in these expressions are well-defined for er = Re(s) > 0, because |e _2s ( r * _r *)| = e~ 2a ^ r *~ r ^ < 1 when r* > r* and 
because the potential Vt(r) decays at least as fast as 1 /r 2 for r —► 00. With these observations in mind it is not 


(27) 

(28) 





difficult to verify that the sequences T± s l(r) obtained by applying k times the operators T± s to the constant function 
£ = 1, converge for all Re(s) > 0 with s / 0 and all r > r c , uniformly on compact intervals, and that ip±(s, •) are 
solutions of Eq. pTj ) fulfilling the required boundary conditions (25). Furthermore, the functions ip±( s i r ) are analytic 
in s for any fixed r > r c . For more details on these assertions we refer the reader to Ref. [33] or Sec. XI .8 in Ref. J3J. 

Next, let us discuss the analytic continuation of the function ip + (s,r) for Re(s) < 0. In this case, the integral 
in Eq. (27) does not converge anymore unless Ve(r) decays exponentially fast. However, if the effective potential Vi 
and the function Af in the definition of the tortoise coordinate r* possess appropriate analytic continuations on the 
complex r plane, ip+(s,r) can be analytically continued to Re(s) < 0 by deforming the path of integration in the 
definition of in Eq. (27). The basic idea, which has been used in Ref. 35] hi the context of the Regge-Wheeler 


equation, relies on the following observation: for s = \s\e lv> and r( — r* = pe la , p > 0, we still have \e~ 2sl ' r *~ r -> \ < 1 
as long as Re(s(r* — r*)) = |s|pcos(<^ + a) > 0. For the integration path in Eq. (27), r( and r* are real and r' > r* 
and consequently, a = 0 which implies that only those frequencies s = pe lv lying in the range \ip\ < 7t/ 2 (that is, 
Re(s) > 0) are admissible. However, choosing a new integration path such that a = —7t/ 2 leads to the admissible 
range 0 < tp < tt, so that the integral in Eq. (27) converges for all Im(s) > 0 provided the analytic continuation of Vn 
decays fast enough along the path (decay equal to or faster than l/|r| 2 is enough). Due to Cauchy’s integral theorem 
the new integration path does yield the same value for (T+ S £)(r) as the one computed using the original path in the 
intersection of the two domains Re(s) > 0 and Im(s) > 0, so by deforming the path in this way we obtain the required 
analytic continuation of i /q_(s,r) on the upper half plane Im(s) > 0. 1 

In our calculations, we choose the following integration path for T +s : 

70(A) = r + e la \, A > 0 , 
with angle a slightly larger than —7r/2, and set 


(T + ,OH = 1 + lj 


1 — exp — 2 s 


dr" 


Af(r") 


Ve(r')£(r')dr ', Re(r) > r c 


(29) 


where it is understood that the integral from r to r' in the exponential is performed along the path 7 a . The analytic 
continuations of the functions Af and V) to complex r and their properties are discussed in the appendix. For large 
|r| and Re(r) > r c , Ve decays at least as fast as l/|r | 2 and Af converges to a positive real constant, so that — r' is 
approximately proportional to r' — r for large |r'|. Hence the integral converges for all Im(s) > 0, as explained above. 

The analytic continuation of the function for Re(s) < 0 can be obtained using similar ideas, 


(T-.OM = 1 -If 


1 — exp 2 s 


dr' 1 


Af(r") 


Ve(r')Ur')dr ', 


(30) 


with r an integration path connecting r with r c . However, in this case particular care has to be taken regarding the 
relation between the tortoise coordinate and the physical radius close to the sonic horizon r = r c , where the function 
1 /Af has a pole. In order to motivate our choice for the integration path T, we approximate 


1 


1 


1 


Af(r) Af'(r c ) r — r c 

for r close to r c . Note that Af'(r c ) > 0 is positive since the surface gravity associated with the sonic horizon is 
positive. As a consequence of the residual theorem, the integral over 1/A f increases by a factor of 2'Ki/Af' (r c ) after 
each revolution along a closed path that winds counter-clockwise around r = r c . In the exponential in the integrand 
on the right-hand side of Eq. (30) this would give rise to a multiplicative factor exp( 47 r«s/A/’ , (r c )) 1 which is bounded 
for all Im(s) > 0 . 

Motivated by these observations, we choose the integration path 


F / 3 (A) = r c + (r - r c )exp(-e I, 3 A), A>0 


1 A similar analytic continuation can be obtained on the lower half plane by choosing a = +7r/2. However, because the functions J\f and 
Vt in Eq. ( |21| ) are real the quasi-normal frequencies s come in complex conjugate pairs, and thus it is sufficient to consider the upper 
half plane. 
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with /3 slightly larger than —7 t/ 2, which spirals counter-clockwise around the point r = r c , see Fig. [2] Along this path 
we have, for s = \s\e l(p , 


exp 2s 


dr h 




N{jr") 


~ exp 


2s 


A f'{r c ) 


(—e ,/3 )dA = exp ( —2|s|e' 


dP+v). 


M'(r c ) 


which is bounded provided |/9 + <^| < 7 t/ 2. Therefore, choosing /3 = —k/2 + 5 with small <5 > 0 guarantees convergence 
of the integral in Eq. (301 for all 0 < < 7r — 5, so by choosing S > 0 small enough we can cover the whole upper 


plane Im(s) > 0. More details and rigorous justifications of our method will be provided elsewhere 



FIG. 2: The integration path Fp in the complex r-plane 


For given Im(s) > 0 we numerically compute the functions ip±(s, r) and their first derivative ip±(s , r) by truncating 
the iteration in Eq. (f26j) to some finite k and computing the operators T± s using Eqs. (29) and (30), where we discretize 


the integrals using the trapezoidal rule. We choose a = —1.57 and j3 = —1.5, and we find that in practice only about 
k ~ 10 iterations are required for good accuracy. The functions J\T(r) and V(,[r) are computed from Eqs. ( [22p3| , 
where n(r) is determined numerically by solving Eq. § via a standard Newton algorithm Ezj. In order to find the 
quasi-nornral frequencies we match the two solutions y>+ and ip- by finding the zeros of their Wronski determinant, 


W(s) := det 


( ip+{s,r ) ip-(s,r) 
\Afip' + (s,r ) A/"*//_ (s,r) 


^ = £+(s,r)A/'£_(s,r) - A7£+(s, r)£_(s, r) + 2s£+(s,r)£_(s,r), 


(31) 


where £±(s,r) = limfc_ >00 (T| s l)(r), at some intermediate point r m (r c < r m < oo) which we typically choose to be 
about r m ~ 1.5r c . The zeros of W are obtained numerically using a standard Newton algorithm m , where the 
derivative of W ( s ) with respect to s is approximated using a simple finite difference operator. 

To test our algorithm, we have applied it to the computation of the quasi-normal frequencies for odd-parity linearized 
gravitational perturbations of a Sclrwarzschild black hole, in which case the functions A f and Vt in Eq. (21) are replaced 
by A f(r) = 1 — rg/r and Vy(r) = £(£ + 1 )/r 2 — 3r#/r 3 , respectively. In the quadrupolar case we found the following 
frequencies: s ■ r H = -0.17792 + 0.747341, -0.54783 + 0.693421, -0.95656 + 0.60211*, -1.4103 + 0.50301*, -1.8937 + 
0.41503*, —2.3912 + 0.33859*, —2.8958 + 0.26651*, which agree to high accuracy with those obtained from Leaver’s 
continued fraction method [221 ■ In order to produce these results we have chosen r m = 1.5r#, discretized the integrals 


in Eqs. ( 29|30 l using 40,000 points and performed 14 Banach iterations. We have varied these numbers in order to 
obtain five significant figures in all the frequencies. 


IV. QUASI-NORMAL OSCILLATIONS FROM A CAUCHY EVOLUTION 


In this section, we solve the Cauchy problem for the wave equation (14) numerically, starting with a Gaussian 
pulse with zero velocity as initial data. We show that a static observer, after registering an initial burst of radiation, 
measures a ringdown signal whose frequency is given by the one of the fundamental quasi-normal mode. 
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A. Reduction to a first-order symmetric hyperbolic system 


We formulate the Cauchy problem for Eq. (141 on the t = const, hypersurfaces of the metric (16) outside the sonic 
horizon. To this purpose we first write the acoustic metric in its ADM form 

© = — a(r) 2 dt 2 + y(r) 2 (dr + f3(r)dt) 2 + R(r ) 2 (dd 2 + sin 2 ddip 2 ) , (32) 

with the functions a(r), /3(r) and y(r) given by 

Ny/N + (u r ) 2 \u r \ 


in N 
&(r) = \ T v s—fr7 > 


h s Vy 


P(r) = 


Y 


.nlVY 
= -Jr- 


where Y := N + (1 — v 2 )(u r ) 2 and \u r \ = |/i|/(r 2 n). Using the following decomposition of 4/ into spherical harmonics, 

r z ' 

lm 

and introducing the auxiliary fields (suppressing the indices lm in what follows) 


7T := - ( d t (j) - Pd r (/)) , 
a 


x ■■= d r <t> , 
7 


Eq. (14) can be cast into first-order symmetric hyperbolic form: 

d t( j) = acfr + 'yPx, 

dtX = ~d r (a Tr + 'yPx),, 

7 


with the effective potential Ue(r) given by 


/ r \ 2 _ 

fR\ 2 , n x 

(s) d - 

( — 1 (olx + 7P7T) 


- aUe(r)q i, 


(33) 

(34) 

(35) 


Ui(r) = 


ajR 2 


a'y 


+ 


1(1 + 1 ) 

R 2 ' 


Explicit evaluation of this potential leads to 

2 


Ui(r) = m 2 


(3n 2 -l)£;(r) + 2K) 2 1 + --(1 -v 2 + W) 

2 n 


i(i+ 1 ) 

R 2 ’ 


(36) 


(37) 


where we recall that E := ?'jy/(4r) — ( u r ) 2 and u r = fj,/(r 2 n). As before, n!/n can be computed using Eq. (24) for all 
r r c and at r = r c we can use th e express ion in Eq. (13) instead. 

We solve the first-order system ( 33|34|35 ) using a finite-difference code based on the method of lines. The spatial 
domain is a finite interval r £ [r c ,r out J with r out r c large enough such that spurious reflections from the outer 
boundary do not affect the wave signal measured by the static observer for the times used in our simulations. There 
are no boundary conditions that must be specified at the inner boundary r = r c since there all the characteristic 
velocities 


ot N 

A 0 = 0, A± = 0 ± - = - 

7 y 


±v s N - (1 - v 2 )y / N+ (u r ) 2 u r 


are zero or positive. At the outer boundary r = r ou t there is one incoming mode 

v in = R(n + x) 

which we set to zero. Although this boundary condition is not exactly transparent to the physical problem, it yields 
only small spurious reflections when r ou t ru and as mentioned above, we extract the physical information only at 
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events which are causally disconnected from the boundary surface in order to make sure that there is no influence 
from the boundary. 

The spatial operators d r are discretized using a fifth-order accurate finite difference operator D 6 _ 5 satisfying the 
summation by parts property and the no-incoming boundary condition is implemented through a penalty method. 
The time derivatives <9 t are discretized using a standard fourth-order Runge-Kutta algorithm. For more details on the 
definition of the -D6-5 operator, the penalty method and numerical time-integrators we refer the reader to Ref. [ 38] 
and references therein. 

We have tested our code for the Regge-Wheeler equation on a Schwarzsc hild back ground metric in ingoing 
Eddington-Finkelstein coordinates [32], for which R = r, y (r) = l/a(r) = a/I + rjj/r, /3(r) = ny/(r 7 (r) 2 ) and 
Ue(r) = —3rH/r 3 + £(£+l)/r 2 are substituted into Eqs. ( 33|34|35 ). We checked fifth-order self-convergence of the field 
cA , and by measuring the wave forms seen by a static observer at r = 20 tr we reproduced the following quasi-normal 
frequencies: s-tr = —0.178 + 0.747* for i = 2, s tr = —0.18541 + 1.19889* for i = 3, and s tr — —0.1883 + 1.61836* 
for t = 4, which agree with those given in the literature, see for example table 2 in Ref. m- 


B. Wave forms for a static observer 


In Fig.[3]we show the time evolution of the acoustic perturbations measured by a static observer located at r = 50 tr 
outside the sonic horizon at r c = 7vr. The initial data for the evolution consists of a Gaussian pulse with zero initial 
velocity, 


/(r) = A exp 


1 tr - ro 

2 \ w 


2 


(38) 


with amplitude A = 1.5, width w = 5.0 tr and centered at tq = 15 tr, and 

0(0, r) = f(r), x(0 ,r) = —/'(r), 7 r( 0 ,r) = 

7 (r) 


J3(r) 

a(r) 


/'(+)■ 


(39) 


In our simulations, we placed the outer boundary at r = 1300*+* and used 2 k x 4000 grid points, where we varied 
k over 0,1,2, 3,4 in order to perform convergence tests. We used a Courant factor of 0.5. The background fluid 
describing the Michel flow is a polytrope with adiabatic index 7 = 1.3333. The quantities shown in the plots of Fig. [ 3 ] 
are the multipolar components of the density contrast , defined as 



where we use Eq. (33) and the definition of the auxiliary field x in order to rewrite dt4>im = oar+/3'yx an d dr<t>tm = 7X, 
respectively. As is apparent from these plots, there is an initial burst of radiation which is followed by several cycles of 
oscillations. The plots corresponding to the cases € = 0,1,2 show that these oscillations are taken over by a power-law 
decay at late times. For the remaining cases i > 2 this is probably also true; however, obtaining the power-law tail 
would require much higher resolution in this case. 

For each t > 0, there is a clear ringdown signal, and we determined the frequency and decay rate of the corresponding 
fundamental quasi-normal oscillations by fitting the numerical data to the function 


Ce at sin (tat — S) 


with free parameters C, <7, w and S. The fit is performed in a time window [£1,^2] where the quasi-normal ringing is 
apparent. The resulting frequencies s = a + ioj are shown in table [I] Since there is no clear ringdown signal for the 
particular case l = 0, only results for i > 0 are shown. The number of significant figures shown has been estimated 
by varying the time window and by comparing the results from different resolutions. 


£= 1 

£ = 2 

1 = 3 

£ = 4 

£= 5 

—0.0080 + 0.0170* 

-0.0081 + 0.0301* 

-0.0081 + 0.04271* 

-0.0081 + 0.0552* 

-0.00811 + 0.0677* 


TABLE I: The quasi-normal fundamental frequencies for r c = 7 th and ^=1,2, 3,4, 5 obtained from the data shown in Fig. [ 3 ] 

As mentioned above, the late time behavior is characterized by a power-law decay, ry ~ t ~ p , as is apparent from 
the plots in Fig. [3] for l = 0,1,2. We have determined the power p, again using a standard fitting routine, obtaining 
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Density contrast for 1=0 



t/r H 

Density contrast for 1=2 



t/r H 

Density contrast for 1=4 



t/r H 


Density contrast for 1=1 



t/r H 

Density contrast for 1=3 



t/r H 

Density contrast for 1=5 



FIG. 3: The density contrast parameter p vs. time t measured by a static observer at r = 50r# for different values of the 
angular momentum parameter £. In all plots, the sonic horizon is located at r c = 7th- Note that only a few oscillations appear 
in the monopolar case £ = 0, preventing us from reading off the quasi-normal frequency in this case. For the cases < = 0,1,2 a 
late time power-law decay is also visible. 


the following results 2 : p = 3.89 ± 0.03 for i = 0, p = 6.62 ± 0.37 for £ = 1, and p = 9.16 ± 0.24 for 1 = 2. The error 
has been estimated by performing the fit in different time windows lying between 3000ryy and 8000r# and by using 


2 See Ref. [40] for a general discussion on the late-time tail decay for wave propagation on a curved spacetime, where it is shown that for 
a certain class of problems the decay only depends on the asymptotic properties of the effective potential. Since in our case the effective 
potential NVi as a function of the tortoise coordinate r* decays as [1(1 + ])/V 2 + Clog(r«)/r®] for large r* with C a nonvanishing 
constant, it follows from the results in Ref. E23 that the fluid potential ^ should decay as t ( 2 ^+ 3 ) for £ > 1. We have verified that 
the function ^ in our simulations reproduce this decay rate for t — 0,1,2 to high accuracy. However, the results in sa do not apply 
directly to the density contrast, which is a nontrivial linear combination of \I/ and its first derivatives. 















13 


different resolutions. 



0 1000 2000 3000 4000 5000 6000 7000 8000 

t/r u 


FIG. 4: Self-convergence test for the case t = 2 in Fig. [3] The top curve corresponds to the error between the results using 
2° x 4000 and 2 1 x 4000 grid points, the second curve to the error using 2 1 x 4000 and 2 2 x 4000 points, etc. The convergence 
factor has been estimated to lie close to 5, indicating fifth-order self-convergence. 


In order to check the validity of our numerical results we have performed several self-convergence tests. In Fig. [4] we 
show a particular example in which we plot the difference of the density contrast function rie m between two consecutive 
resolutions. This plot corresponds to the quadrupolar case £ = 2 in Fig. [3] Note that there is high frequency noise 
appearing at £ ~ 1500rqf, which is probably due to the presence of the time and space derivatives of 4> in the expression 
for r]e m in Eq. (40). However, it is clear from the plot that the error decreases with increasing resolution. We have 
estimated the convergence factor to lie close to five, indicating fifth-order self-convergence. 


V. RESULTS FOR THE QUASI-NORMAL FREQUENCIES 


In this section, we present and analyze the results from our calculations of the quasi-normal acoustic frequencies as 
a function of the sonic radius r c and the angular momentum £. All the calculations in this section refer to the Michel 
flow on a Schwarzschild background for a polytropic fluid with adiabatic index 7 = 1.3333. In Sec. |V A[ w e disc uss 
the fundamental frequencies for values of r c ranging in the interval [2r#,30rff] and £ = 0 , 1 , 
also discuss quasi-normal frequencies corresponding to the first few overtones. 


,7. In Sec. VB 


we 


A. Fundamental frequencies 


In table [TT] we show the fundamental monopolar, dipolar and quadrupolar quasi-normal frequencies for different 
values of r c . These frequencies were calculated using the matching method described in Sec. |III[ and in the dipolar and 
quadrupolar cases with r c /r# = 2, 7,10, 20, 30 also using the numerical Cauchy evolution described in the previous 
section. As can be seen from this table, the two approaches give results which are consistent within their numerical 
errors. 

Also shown in table [II] are the values for the surface gravity k of the acoustic metric, computed using Eq. (20 1 . It 
turns out k plays an important role for understanding the behavior of the quasi-normal frequencies as a function of 
the location of the sonic horizon r c . Indeed, k has units of frequency (in geometrized units) and thus it is natural 
to analyze the quasi-normal frequencies in units of k . In Fig. [ 5 ] we show plots of s/k vs . r c for the fundamental 
quasi-normal frequencies s = a + iuj for different values of l. As is apparent from these plots, the value of s/k seems 
to be almost independent of r c for large r c /rn- Specifically, we have found that the empiric formula 


s 

K 


10 < — <30, £=1,2,..., 7, 

rH 


0.387 + (0.21 + 0.606£)i, 


(41) 
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r c /r H 

K • r H 

s ■ r H (£ = 0) 

s-r H (£ = 1) 

s-r H (£ = 2) 

2 

0.16536 

—0.05947 + 0.02661* 

-0.06174 + 0.1398* 
(-0.06 + 0.14*) 

-0.06203 + 0.2416* 
(-0.062 + 0.242*) 

3 

0.08334 

-0.02932 + 0.008119* 

-0.03144 + 0.06813* 

-0.03162 + 0.1191* 

4 

0.05182 

-0.01805 + 0.003508* 

-0.01965 + 0.04204* 

-0.01977 + 0.07381* 

5 

0.03606 

-0.01249 + 0.001812* 

-0.01372 + 0.02920* 

-0.01380 + 0.05138* 

6 

0.02690 

-0.009286 + 0.001042* 

-0.01026 + 0.02178* 

-0.01032 + 0.03838* 

7 

0.02104 

-0.007250 + 0.0006431* 

-0.008036 + 0.01704* 
(-0.0080 + 0.0170*) 

-0.008086 + 0.03006* 
(-0.0081 + 0.0301*) 

8 

0.01703 

-0.005858 + 0.0004166* 

-0.006513 + 0.01381* 

-0.006553 + 0.02436* 

9 

0.01414 

-0.004861 + 0.0002792* 

-0.005416 + 0.01148* 

-0.005449 + 0.02026* 

10 

0.01199 

-0.004118 + 0.0001913* 

-0.004595 + 0.009736* 
(-0.0046 + 0.0097*) 

-0.004623 + 0.01720* 
(-0.00462 + 0.0172*) 

20 

0.00410 


-0.001577 + 0.003344* 
(-0.0016 + 0.0033*) 

-0.001586 + 0.005914* 
(-0.0016 + 0.00591*) 

30 

0.00220 


-0.0008493 + 0.001803* 
(-0.00085 + 0.0018*) 

-0.0008546 + 0.003190* 
(-0.00085 + 0.00319*) 


TABLE II: Fundamental quasi-normal frequencies for acoustic perturbations of the Michel flow for different values of r c and £. 


The frequencies in the first line of each entry for r c /rH are the ones obtained from the matching procedure discussed in Sec. Ill 
and four significant figures are shown. The frequencies in parenthesis refer to the ones obtained from the Cauchy evolution 
code and are shown for comparison. In the monopolar case £ = 0 we have not been able to obtain the frequencies from the 
Cauchy evolution, for the reasons described in the previous section. For £ = 0 and r c > 15rtr we have not been able to compute 
the frequencies in a reliable way using our matching procedure; their computation seems to require higher accuracy than the 
one available in our current code. 


gives a fit for the fundamental frequency to a relative accuracy better than 2%. Notice that in the monopolar case 
l = 0 the behavior of a as a function of r c is different than for higher multipoles t > 1. 




e 


- 0.34 


- 0.35 


- 0.36 


- 0.37 


- 0.38 


- 0.39 



0 5 10 15 20 25 30 

r c /r H 


10 15 20 25 30 

r c /r H 


FIG. 5: The fundamental quasi-normal frequencies in units of k as a function of r c . Left panel: real part a / k divided by 
the surface gravity k for £ = 0,1, 2,..., 7. As is apparent from the plot, for < ^ 0 and r c /rH > 10 these values are almost 
independent of £ and r c , and can be approximated by —0.387 to about 1% accuracy. Right panel: imaginary part cj/k divided 
by the surface gravity for £ = 0,1, 2,..., 7. These values are almost independent of r c , and we found that for £ ^ 0 they are 
well-approximated by the empiric formula 0.21 + 0.606A 


B. Overtones 


Using the matching procedure described in Sec. |III| we have also computed the quasi-normal frequencies of the first 
few excited modes. In table III we present two examples for the quasi-normal spectrum, referring to dipolar and 




















15 


quadrupolar acoustic perturbations, respectively, with r c = 10 rjj- 


n 0 

s • r H (t = 1) 

s ■ r H (i = 2) 

i 

-0.01503 + 0.007955* 

-0.01437 + 0.01590* 

2 

-0.02691 + 0.006537* 

-0.02526 + 0.01408* 

3 

-0.03907 + 0.005732* 

-0.03699 + 0.01257* 

4 

-0.05123 + 0.005219* 

-0.04905 + 0.01149* 

5 

-0.06336 + 0.004855* 

-0.06120 + 0.01071* 

6 


-0.07336 + 0.01012* 

7 


-0.08551 + 0.009659* 

8 


-0.09764 + 0.009281* 

9 


-0.1098 + 0.008965* 


TABLE III: Quasi-normal dipolar and quadrupolar frequency spectrum for acoustic perturbations of the Michel flow with sonic 
horizon located at r c = 10 th- n 0 = 1 denotes the first overtone, n 0 = 2 the second etc. Four significant figures are shown. 
Modes with excitation numbers no > 5 for l = 1 and no > 9 for i = 2 could not be obtained in a reliable way with the current 
version of our code, since the computation of their frequency seems to require a more powerful Newton algorithm or higher 
accuracy. 

In Fig. [6] we show plots of the quasi-normal dipolar and quadrupolar spectrum in units of the surface gravity « for 
different values of r c . As in the case of the fundamental frequencies, we appreciate from these plots that the spectrum 
of quasi-normal excitations seems to be nearly independent of r c , indicating that the frequencies scale like k. 

1.6 

1.4 

1.2 


* 

0.8 

0.6 

0.4 

0.2 

-10 -8 -6 -4 -2 0 

o/k 

FIG. 6: The spectrum of the quasi-normal acoustic excitations. Shown are the imaginary vs. the real part of the frequencies 
s/k, divided by the surface gravity for r c /rjf = 2,5,10 and £=1,2. As is apparent from the plot, the spectrum is approximately 
independent of r c . 


r c /r H = 5,1= 1 
1 -^ = 10 , 1=1 
r c /r H = 2,1=2 
r c /r H = 5,1=2 
1 ^ = 10 , 1 = 2 


C. Eikonal limit 

In the high-frequency limit, the quasi-normal oscillations can be interpreted in terms of wave packets which are 
concentrated along a circular null geodesics and decay because the circular null geodesic is unstable, see Refs. GHEE] 
and references therein for more details. Therefore, one expects that in this limit the quasi-normal frequencies s are 
related to the properties of the unstable circular null geodesics. As shown in [41] the imaginary part u> = Im(s) of 
s, describing the oscillatory behavior, is directly related to the angular velocity of the unstable circular null geodesic, 
while the real part a = Re(s) is equal to its Lyapunov exponent. 

As can be deduced from the analysis in m an arbitrary asymptotically flat, static spherically symmetric metric of 
the form 


ds 2 


-f(r)dT 2 


+ r 2 (dd 2 + sin 2 ddip 2 ) 

9{r) 


(42) 
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with time coordinate T and positive smooth functions f(r) and g{r) possesses an unstable circular null geodesic at 
r = v c irc if and only if the function H(r ) := f(r)/r 2 has a local maximum at r = r c i rc , and in this case the associated 
angular velocity and Lyapunov exponent are given by 




— \/ H ( t'cir'c )' A circ — 


f(r)g(r) 


1 d 2 
H(r ) dr 2 


H (r) 


(43) 


For high values of £, these parameters determine the quasi-normal frequencies according to the formula 


s — (n 0 l/2)X C i rc i££l c 


(44) 


with n a the overtone number, see HU- Comparing Eq. ( |'12[ ) with the form ( pL9[ ) of the acoustic metric and discarding 
the conformal factor n/(hv s ) which does not affect the null geodesics as trajectories in spacetime, we find that in our 
case /(r) = X(r)v 2 and g(r) = X(r), such that 


X(r)v 2 

H{r) = 


VK r )g( r ) = x{r)v s . 

A plot of the function H(r ) for the case r c /rH = 10 is given in Fig. [7j which shows that the acoustic metric 



FIG. 7: Graph of the function H(r) for the case where the sonic horizon is located at r c /rjj = 10. As is clearly visible from 
the plot, this function has a maximum where the acoustic metric has an unstable circular null geodesics. This maximum is 
numerically determined to be located at r c i re /rir = 14.158. 


for the Michel flow admits unstable circular null geodesics. Numerically, we find the values r C i rc /rH = 14.158, 
^circ/ft = 0.59159 and A c j rc /(2K) = 0.38637 which agree remarkably well with the corresponding values in the 
empirical formula (411 describing the fundamental frequencies. We have repeated the analysis for higher values of 
r c /ru ranging between 10 and 30, finding similar values for f l C ir C /^ and A c i rc /(2K) (the difference is less than 1%). 


VI. CONCLUSIONS 

In this work, we have analyzed spherical and nonspherical acoustic perturbations of the Michel flow, which describes 
a perfect fluid which falls radially into a Schwarzschild black hole. As shown by Moncrief m , the equations of motion 
for such perturbations can be cast into a wave equation on a curved effective background geometry described by the 
acoustic metric. For the case of the Michel flow, the acoustic metric has the same qualitative properties as a black 
hole spacetime and thus describes a natural analogue black hole. 

Using this natural astrophysical analogue black hole, we have shown by numerical computation that when perturbed, 
the Michel flow exhibits quasi-normal acoustic oscillations. We have computed the associated frequencies s = a + ioj 
using two different methods. The first method which, to our knowledge, is new is based on matching the two local 
solutions ij}+(s, r ) and ip~{s, r) of the radial mode equation which, for Re(s) > 0, are decaying as r —> oo and r —> r c , 
respectively. A common challenge for computing the quasi-normal modes is to determine the analytic continuation 
of these functions for Re(s) < 0 and to find the complex frequencies s for which ip + and if}- are linearly dependent. 
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While in some cases the solutions ip± can be represented by simple series expansions and the quasi-normal frequencies 
can be found using continued fraction techniques |22j . in our problem the effective potential appearing in the mode 
equation is not even known in closed form and so more general methods are required. The new ingredient of our 
method consists in computing the analytic continuations of ip± via a Banach iteration technique, where each iteration 
leads to an improved approximation for the solution. Each iteration involves computing a line integral in the complex 
r-plane which converges for all to = Im(s) > 0. While the integral in each iteration needs to be computed accurately, 
we have found that only a few iterations are needed in order to achieve high accuracy. The two solutions are then 
matched by finding the zeros of their Wronski determinant using a standard Newton algorithm. 

Our method is rather general and does not depend on the details of the effective potential except for the fact that it 
should possess a sufficiently well-behaved analytic continuation on the complex /'-plane. What precisely we mean by 
“sufficiently well-behaved” will be explained in detail elsewhere [3B], but it seems flexible enough to comprise many 
relevant effective potentials found in general relativity (including the Regge-Wheeler potential and its generalization to 
the Reissner-Nordstrom case). Another advantage of our method is that it does not require a closed-form expression 
for the effective potential. For the case of acoustic perturbations of the Michel flow considered in this article the 
potential is only known in implicit form, though it is analytic in 1/r as we have shown in the appendix. In the cases 
we have analyzed here, our method seems to work very well to find the fundamental frequencies and the first few 
overtones. However, so far our code fails to find very high overtones. The reason for this is probably related to our 
simple Newton algorithm and our crude finite-difference approximation for the derivative of Wronski determinant. 

Using our method we have computed the quasi-normal acoustic frequencies of the Michel flow for different values r c 
and ru of the sonic and event horizon radii, and for different values of the angular momentum number £. By means 
of the Cauchy code described in Sec. [IV]we have verified the validity of the fundamental frequency for £ > 0, and also 
computed the late time power-law decay rate in some cases. Although in general the frequency spectrum depends 
on two parameters r# and r c , or, equivalently, on rn and the surface gravity k of the acoustic hole, we found that 
for r c th the quasi-normal frequencies s scale like «, the parameter th becoming unimportant. Furthermore, for 
r c rn the real part of s describing the decay rate depends only mildly on £ for £ > 1. Specifically, we have found 
the following empiric formula for the fundamental frequency: 

sk- 1 ~ -0.387+ (0.21 + 0.606A)* (45) 


for £ = 1, 2,..., 7 and r c ranging in the interval between 10 th and 30ry/, with rjj the event horizon radius. In the 
limit where the sound speed Voo <C c at infinity is much smaller than the speed of light, n can be given by a simple 
analytic formula. It follows from Eq. (20) and standard expansions in v ^ := v <*, / c mm that 



for a polytropic equation of state with 7 = 4/3. For a black hole of mass M this gives 

.3 


k~8x!0 5 


Mq 

M 


OO 

s 


(46) 


(47) 


with Mq the solar mass. 

Although in this work we have restricted ourselves to a polytropic fluid with adiabatic index 7 = 1.3333 ~ 4/3, 
other fluid flows could be analyzed with our method, provided they are described by an analytic equation of state 
satisfying the assumptions (F1)-(F3) listed in Sec.[TTj Furthermore, based on our general results in Ref. HU, it should 
not be difficult to generalize our calculations to more general nonrotating black holes, and to analyze the dependency 
of the quasi-normal acoustic frequencies on the background metric. It would be interesting to study the impact of 
these acoustic oscillations on the emission of electromagnetic and gravitational radiation. 
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Appendix A: Analytic continuation of the functions N and Vi 


In this appendix, we prove that the functions Af{r) and V( (r) in the mode equation (21) admit analytic continuations 
on the domain Re(r) > r# with the properties that 


lim A f(r) = v a 

1 ->-00 V 7 

Re(r)>r H 


r iiin r 2 V e (r) = v 00 i{l + 1), 
Re(r)>rjr 


(Al) 


where := v s (noo) > 0 is the sound speed at infinity. For this, we need to assume that in addition to the properties 
(F1)-(F3) the specific enthalpy h(n) is an analytic function of n. For definiteness, we shall assume that h(n) is given 
by the polytropic equation of state, Eq. which is analytic on the domain Re(n) > 0. 

Under these assumptions, we first prove that the Michel flow solution n(r), which is implicitly determined by 
Eq. possesses an analytic continuation on the domain Re(r) > m such that 


lim n(r) = n a 
1 ->-00 V 7 

Re(r)>r/f 


(A2) 


where n oo > 0 is the particle density at infinity. In order to prove this statement, following m we introduce 
dimensionless quantities x := r/rjj , z '■= n/no, no := (eo/AT) 1 ^ 7-1 ', in terms of which Eq. (I6j) can be rewritten as 


F„(x,z) := f(z)- 


1-i 


F 


= foo= const., 


(A3) 


where f(z ) = 1 + z 7 1 = l + e^ 7 1 U°g( z ) j s the dimensionless enthalpy function and /oo = /(z^), -oo > 0, its value at 


infinity. The function F ^ defined by Eq. (A3) is analytic on the domain f l c := {(x,z) £ C 2 : Re(x) > 0,Re(z) > 0}. 


In m we showed that there exists a unique real-valued differentiable function zq : [1, oo) — > R. (the Michel solution), 
defined on and outside the event horizon, such that F ll (x,zo(x)) = for all x > 1 and lim x _>.oo z{x) = z^. This 
solution has the property that the partial derivative of F,, with respect to z, 


dF 2 /(z) 2 2 

~^{x,z) = - v{z) 

oz z 


l- 1 —' 1 

X 




-1 


F 


2 1 


(A4) 


is different from zero for all x > 1 except at the location of the critical point x = x c . By continuity, dF^/dz is also 
different from zero in an open neighborhood U C S2 C of the graph G := {(x,zo(x)) : x > l,x ^ x c }. Therefore, it 
follows from the implicit function theorem that for an open neighborhood V C U of G in S2 C , zq{x) admits a unique 
analytic continuation z(x) whose graph lies in V and such that F fl (x,z(x)) = /£, for all (x,z(x)) £ V. 

It remains to prove that z(x) can be further extended to a neighborhood of x = oo and to an open neighborhood 
of the critical point. For the former case, we introduce the new variable y := 1/x and rewrite Eq. (A3) as 


Fn(y,z) : = f h\ 7> z ) = /(A) 


1 , y 2 4 

1 - 2 / + -tv 

Z z 


= f 2 

J o 


The function F ^ is analytic on the domain y G C, R e(z) > 0, and it satisfies F /i (0, Zqq) = /£> and 

dF 2fl 2 

o (0? -^oo) — ^(^oo) T 1 0. 

UZ Zqq 

Therefore, it follows from the implicit function theorem that there exists an open neighborhood V of (0, z^) and 
a unique function z(y) whose graph lies in V such that 5(0) = z^ and F^(y,z(y)) = /£, for all ( y,z(y )) £ V. By 
uniqueness of the analytic continuation, z(x) = z( 1/x) for large enough \x\, which proves that the analytic extension 
of z(x) exists for sufficiently large \x\. Furthermore, 


lim z(x) = lim z(y) = Zoo- 

x^-oo v ' v ^o 

Re(®)>o 

Next, we discuss the analytic continuation of Zq(x ) in an open neighborhood of the critical point x = x c . For this, 
we first note that in a vicinity of the critical point (x c ,z c = zo(x c )) the function Fy has the Taylor representation 


F^ (x c + £, z c + C) — Fy [x c , z c ) T 


&F Uf _ _ ,, 2 , _ ,,, , d 2 F 

[x c , z c )£ +2 (x c , z c )/C + „ (a 


dx 2 


dxdz 


dz 2 


c)C 2 


+ ^3(C) 0, 
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where the error term i? 3 (£, () is at least cubic in (£,£). Let z' c € R denote one of the two roots of the quadratic 
polynomial (cf. Eq. ©) 


d 2 F M 
dx 2 


g 2 p 


g 2 f 

c )4 + ^(xc 


^c)(4) 2 = 0, 


and introduce the function 




j 2 [i^c + + z'^rf) - F^(x c , z c )] 


dx a* (^c) ^c) + 2 (x c , z c )z c r] ■ 


9 2 G 




for £ ^ 0 , 
for £ = 0 . 


Then, H „ is analytic in an open neighborhood of (£, 77 ) = (0,1) in C 2 , satisfies iJ M (0,1) = 0 and 




dr] 


dxdz 


B 2 F 

z c) z 'c + ~Q^r( x c, z c ){z' c ) 2 = =f3 


h 2 c y/1 + 3{yl - W e ) 
xl2± v/1 + 3(i/2 - W c ) 


7 ^ 0. 


Therefore, using once again the implicit function theorem, it follows the existence of an open neighborhood Z of (0,1) 
in C 2 and a unique analytic function r]{£) whose graph lies inside Z such that 77 ( 0 ) = 1 and t?(£)) = 0 for all 

(Ci ??(£)) £ Z. By construction z(x) := z c + z' c (x — x c )r](x — x c ) is analytic and satisfies F^(x, z(x)) = F(x c , z c ) = 

This demonstrates the existence of the analytic continuation of z(x) in a neighborhood of the critical point. 

With these results, it follows directly from Eqs. ( 22|23|24 ) that the functions J\f(r) and Ve(r) have analytic contin¬ 
uations for complex r, and that these continuations satisfy Eq. (Al). 


[1] Event horizon telescope, http://www.eventhorizontelescope.org. 

[2] S. Doeleman and et al. Event-horizon-scale structure in the supermassive black hole candidate at the galactic centre. 
Nature, 455:78, 2008. 

[3] A.E. Broderick, T. Johannsen, A. Loeb, and D. Psaltis. Testing the no-hair theorem with event horizon telescope obser¬ 
vations of Sagittarius A*. Astrophys. J., 784:7, 2014. 

[4] O. Zanotti, J.A. Font, L. Rezzolla, and P.J. Montero. Dynamics of oscillating relativistic tori around Kerr black holes. 
Mon.Not.Roy.Astron.Soc. , 356:1371-1382, 2005. 

[5] M. Megevand, M. Anderson, J. Frank, E.W. Hirschmann, L. Lehner, S.L. Liebling, P.M. Motl, and D. Neilsen. Perturbed 
disks get shocked, binary black hole merger effects on accretion disks. Phys.Rev., D80:024012, 2009. 

[6] M. Anderson, L. Lehner, M. Megevand, and D. Neilsen. Post-merger electromagnetic emissions from disks perturbed by 
binary black holes. Phys.Rev., D81:044004, 2010. 

[7] F.C. Michel. Accretion of matter by condensed objects. Astrophysics and Space Science, 15:153-160, 1972. 

[8] H. Bondi. On spherically symmetrical accretion. Monthly Notices Roy Astronom. Soc., 112:195-204, 1952. 

[9] S.L. Shapiro and S.A. Teukolsky. Black Holes, White Dwarfs, and Neutron Stars. John Wiley & Sons, New York, 1983. 

[10] F.S. Guzman and F.D. Lora-Clavijo. Exploring the effects of pressure on the radial accretion of dark matter by a 
Schwarzschild supermassive black hole. Mon. Not. R. Astron. Soc., 415:225-234, 2011. 

[11] E. Chaverra and O. Sarbach. Radial accretion flows on static, spherically symmetric black holes. 2015. arXiv:1501.01641. 

[12] V. Moncrief. Stability of stationary, spherical accretion onto a Schwarzschild black hole. Astrophys. J., 235:1038-1046, 
1980. 

[13] N. Bilic. Relativistic acoustic geometry. Class. Quantum Grav., 16(12):3953-3964, 1999. 

[14] C. Barcelo, S. Liberati, and M. Visser. Analogue gravity. Living Rev.Rel, 8:12, 2005. 

[151 P. Mach and E. Malec. Stability of relativistic Bondi accretion in Schwarzschild-(anti-)de Sitter spacetimes. Phys.Rev. D, 
D88:084055, 2013. 

[16] T.K. Das, N. Bilic, and S. Dasgupta. Black-hole accretion disc as an analogue gravity model. JCAP, 0706:009, 2007. 

[17] T.K. Das. Astrophysical accretion as an analogue gravity phenomena. 2007. arXiv:0704.3618. 

[18] D.B. Ananda, S. Bhattacharya, and T.K. Das. Acoustic geometry through perturbation of mass accretion rate I - radial 
flow in general static spacetime. 2014. arXiv: 1406.4262. 

[19] D.B. Ananda, S. Bhattacharya, and T.K. Das. Acoustic geometry through perturbation of mass accretion rate - axisym- 
metric flow in static spacetimes. 2014. arXiv: 1407.2268. 

[20] E. Berti, V. Cardoso, and J. Lemos. Quasinormal modes and classical wave propagation in analogue black holes. Phys. 
Rev. D, 70:124006, 2004. 

[21] S.R. Dolan, L.A. Oliveira, and Luis L.C.B. Crispino. Quasinormal modes and regge poles of the canonical acoustic hole. 
Phys.Rev., D82:084037, 2010. 

[22] E. W. Leaver. An analytic representation for the quasi-normal modes of Kerr black holes. Proc. R. Soc Lond. A, 402:285-29, 
1985. 



















20 


[23] N. Andersson and C.J. Howls. The asymptotic quasinormal mode spectrum of nonrotating black holes. Class. Quantum 
Grav., 21:1623-1642, 2004. 

[24] A. Zimmerman, H. Yang, Z. Mark, Y. Chen, and L. Lehner. Quasinormal modes beyond Kerr. Astrophys. Space Sci. 
Proc., 40:217, 2015. 

[25] E. Chaverra and O. Sarbach. Polytropic spherical accretion flows on Schwarzschild black holes. AIP Conf.Proc., 1473:54-58, 
2012 . 

[26] J. Karkowski and E. Malec. Bondi accretion onto cosmological black holes. Phys.Rev. D87 (2013), D87:044007, 2013. 

[27] M. Heusler. Black Hole Uniqueness Theorems. Cambridge University Press, Cambridge, England, 1996. 

[28] R.M. Wald. General Relativity. The University of Chicago Press, Chicago, London, 1984. 

[29] H.-P. Nollert and B.G. Schmidt. Quasinormal modes of Schwarzschild black holes: Defined and calculated via laplace 
transformation. Phys. Rev. D, 45:2617-2627, 1992. 

[30] H.-P. Nollert. TOPICAL REVIEW: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars. Class. 
Quantum Grav., 16:R159-R216, 1999. 

[31] K.D. Kokkotas and B.G. Schmidt. Quasi-normal modes of stars and black holes. Living Reviews in Relativity, 2(2), 1999. 

[32] E. Berti, V. Cardoso, and A.O. Starinets. Quasinormal modes of black holes and black branes. Class. Quantum Grav., 
26:163001, 2009. 

[33] R.G. Newton. Analytic properties of radial wave functions. J. Math. Phys., 1:319-347, 1960. 

[34] M. Reed and B. Simon. Methods of Modern Mathematical Physics, Vol. Ill: Scattering Theory. Academic Press, San Diego, 
1979. 

[35] B.P. Jensen and P. Candelas. The Schwarzschild radial functions. Phys. Rev. D, 33:1590-1595, 1986. 

[36] E. Chaverra and O. Sarbach. In preparation, 2015. 

[37] W. H. Press, S. A. Teukolsky, W. T. Watterling, and B. P. Flannery. Numerical Recipes in Fortran. Cambridge University 
Press, Cambridge, 1992. 

[38] O. Sarbach and M. Tiglio. Continuum and discrete initial-boundary value problems and Einstein’s field equations. Living 
Rev. Relativity, 15, 2012. http://www.livingreviews.org/lrr-2012-9. 

[39] O. Sarbach and M. Tiglio. Gauge invariant perturbations of Schwarzschild black holes in horizon penetrating coordinates. 
Phys. Rev. D, 64:084016(1)-(15), 2001. 

[40] E.S.C. Ching, P.T. Leung, W.M. Suen, and K. Young. Late time tail of wave propagation on curved space-time. Phys. 
Rev. Lett., 74:2414-2417, 1995. 

[41] V. Cardoso, A.S. Miranda, E. Berti, H. Witek, and V. Zanchin. Geodesic stability, Lyapunov exponents and quasinormal 
modes. Phys. Rev. D, 79:064016, 2009. 



